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Irreversible difFusion limited cluster aggregation (DLCA) of hard spheres was simulated using 
Brownian cluster dynamics. Bound spheres were allowed to move freely within a specified range, 
but no bond breaking was allowed. The structure and size distribution of the clusters was investi- 
gated before gelation. The pair correlation function and the static structure factor of the gels were 
determined as a function of the volume fraction and time. Bond flexibility led to local densification 
of the clusters and the gels, with a certain degree of order. At low volume fractions densification of 
the clusters occurred during their growth, but at higher volume fractions it occurred mainly after 
gelation. At very low volume fractions, the large scale structure (fractal dimension), size distribu- 
tion and growth kinetics of the clusters was found to be close to that known for DLCA with rigid 
bonds. Restructuring of the gels continued for long times, indicating that aging processes in systems 
with strong attraction do not necessarily involve bond breaking. The mean square displacement of 
particles in the gels was determined. It is shown to be highly heterogeneous and to increase with 
decreasing volume fraction. 
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I. INTRODUCTION 

Strong attraction between particles in solution leads 
to aggregation. The kinetics of this aggregation pro- 
cess depends on the probability that a bond is formed 
when two particles collide. Two limiting cases are dif- 
fusion limited cluster aggregation (DLCA) if bonds are 
formed at each collision and reaction limited cluster ag- 
gregation (RLCA) if the probability to form a bond is 
very small. Irreversible aggregation has been studied 
in detail both experimentallji^i^i^i^i^ and using com- 
puter simulations'^^iSiiaiiiii^ii^iiiii^i^iil. The structure 
of clusters formed by random aggregation is self simi- 
lar and characterized by a fractal dimension (d/), which 
relates the radius of gyration {Rg) to the aggregation 

number (to): to oc Rg^ ■ The number of clusters with 
aggregation number to can be described by a power law: 
N(rn) oc m~'^ . As long as the clusters are on average 
far apart (flocculation) one finds df = 1.8 and r = for 
DLCA and d/ = 2.1 and r = 1.5 for RLCAi^. 

With time the clusters grow and the cumulated vol- 
ume occupied by the clusters {Vcum = -^(w) •47ri?^/3) 
increases so that the average free space between the ag- 
gregates decreases. When Vcum approaches the volume 
of the system the aggregates start to interpenetrate. The 
aggregation process of highly interpenetrated clusters can 
be described by the percolation model and leads to gela- 
tion. For percolating clusters df = 2.5 and r = 2.2^^. 
The cross-over between flocculation and percolation oc- 
curs at a characteristic aggregation number (toc) and ra- 
dius of gyration (Re), that decrease with increasing par- 
ticle volume fraction (0)^ii£iii^. 

Computer simulations of irreversible DLCA and RLCA 
have been done so far for hard spheres that form rigid 
bonds at contac t^^i^^'^''i^° . In this case the aggregated 
particles are on average bound to two other particles. 



because ternary collisions are not possible. However, in 
reality the bonds may be flexible, i.e. they may freely 
rotate. One example is the much studied aggregation of 
spheres in the presence of other smaller particles through 
a depletion interactio n^^'^^i^^'^'^i^^'^^ . Another example 
is the aggregation of emulsion droplets with a slippery 
layer"'^. The latter experiment has motivated computer 
simulations of irreversible diffusion limited aggregation 
(DLA) with flexible bonds^-. The difference between 
DLA and DLCA is that during the former individual 
particles are allowed to diffuse until they collide with a 
single cluster>i^, while during the latter all particles in the 
systems diffuse and collide to form many clusters. DLA 
leads to a self similar cluster with df = 2.5. In refJ^ 
the DLA simulation was modified to include free diffu- 
sion of the particle on the surface of the particle to which 
it is bound and was called slippery DLA. Real random 
aggregation processes are, of course, better described by 
DLCA. 

Here we report on a simulation study of DLCA with 
finite interaction range in which the relative motion of 
bound particles is unhindered as long they remain within 
each others range. If the interaction range is very small 
compared to the radius of the particles this method could 
be called slippery DLCA, but here we will use the expres- 
sion flexible DLCA for all interaction ranges. Flexible 
DLCA (or RLCA) should represent realistically the ex- 
perimental systems mentioned above in the limiting case 
that the attraction energy is much larger than the ther- 
mal energy. The results of flexible DLCA will be com- 
pared to those obtained by DLCA with rigid bonds. 



II. SIMULATION METHOD 

The simulation method used here is called Brownian 
Cluster Dynamics (BCD). A detailed description of the 
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method and a comparison with molecular dynamics were 
reported elsewhere?^. Briefly, clusters are formed by con- 
necting spheres within each others interaction range with 
probability P. Particles are chosen randomly and moved 
a step size s in a random direction unless it leads to over- 
lap or breaks a bond. The centre of mass displacement 
of the clusters is calculated and the clusters are moved 
cooperatively in the same direction so that the total dis- 
placement is inversely proportional to their radius, unless 
it leads to overlap. BCD is equivalent to molecular dy- 
namics if the cooperative cluster displacement is omitted 
as long as s is sufficiently small. Systems with rigid bonds 
are simulated by performing only the cooperative cluster 
movements and not the individual particle displacements 
within the clusters. The equilibrium state obtained by 
molecular dynamics and BCD is the same, but the dy- 
namics depend strongly on whether the bonds are rigid 
or flexible. 

Irreversible DLCA is simulated by setting P — 1. 
The simulation is started with Ntot randomly distributed 
spheres with unit diameter in a box of size L so that 
(j) — Ntot / (t^ / Q) ■ The unit of time is set equal to the 
time needed for an isolated sphere to diffuse a distance 
equal to its diameter. E.g. for spheres with diameter 
1/Ltm in water at 20°C the time unit is 0.4 seconds. The 
box size was varied up to L = 100 and all the results 
shown here were not influenced by finite size effects un- 
less specified. 



III. RESULTS 

In the following we show mainly results obtained with 
the interaction range fixed at e = 0.1, but we will briefly 
discuss the effect of varying the range. 

A. Kinetics 

It is well known that in the fiocculation regime, 
i.e. Vcum ^ the cluster growth during DLCA can 
be described by the kinetic equations introduced by 
Smolechowski22i2^i2a: 

^lj2K{hj)Nm{j)-Y,K{m,j)N{m)N{j) 

i+j J 

(1) 

Where K{i,j) is the so-called kernel that expresses the 
rate constant at which a cluster with aggregation number 
i collides with a cluster with aggregation number j: 

K{i,j) = in {RcoL, + Rcoi,j) (A + Dj) (2) 

with D the diffusion coefficient of the clusters and Rcoi 
their collision radius. D is inversely proportional to the 
hydrodynamic radius (Rh) of the clusters. For large clus- 
ters both Rcoi and Rh are proportional to Rg^^. In fact, 
the growth kinetics of DLCA can be well described by 



eq([T]) using a constant kernel equal to K{1, 1). It follows 
that the time dependence of the weight (m^,) averaged 
aggregation number can be written as: 

m^, = l + K{l,l)t<j) (3) 

In our simulations the particle diameter is unity and 
the diffusion coefficient of the individual spheres is 1/6 
so that K(l,l) — Att/S. The ratio of the weight and 
the number (m„) averaged aggregation numbers becomes 
two for large m-u, ■ These results have been confirmed us- 
ing BCD with rigid bonds, but only if (jj was very low. 
In^'' the interaction range was zero. Increasing the inter- 
action range increases the collision radius and therefore 
K{1, 1), but the effect is smafl for e = 0.1. For (p > 0.01 it 
is difficult to observe the fiocculation regime because the 
crossover to the percolation regime occurs already when 
is still relatively small. The growth of acceler- 
ates when Vcum approaches L'^ and diverges at the 
gel time (tg). This change of the growth kinetics between 
fiocculation at short times and percolation close to tg has 
also been observed experimentally^. 

FigurelU^) compares the increase of as a function 
of time for DLCA with rigid and with flexible bonds at 
different volume fractions, while fig.((T|D) shows the evo- 
lution of the polydispersity index {raw /ran). The size of 
the clusters at the start of the simulation depends on the 
interaction range and the volume fraction, because more 
particles are within each others range if e or are larger. 
For large e and (p the system percolates immediately, i.e. 
when the spheres are still randomly distributed, as was 
discussed in detail inS. The percolation threshold de- 
creases with increasing (f) and is 0.27 for e — 0.1. 

At (j) = 0.01 the increase of ra^ with time was linear 
and aiw/ran was close to two as expected for the fioccula- 
tion regime. The results were similar for rigid and flexible 
DLCA, but rigid DLCA showed the influence of a tran- 
sition to the percolation regime sooner. The gel point 
was not reached for this volume fraction during the sim- 
ulation which took several weeks of computer time. For 
(j) > 0.08 the fiocculation regime was not observed and 
the growth of both the size and the polydispersity of the 
clusters accelerated until they diverged at tg . The growth 
of the clusters and gelation was somewhat faster for flex- 
ible DLCA. At 4> = 0.05 the increase of ra^ was initially 
slightly faster during flexible DLCA, but then the growth 
became almost linear over a short period of time while 
the polydispersity index remained constant. After this 
period the increase of ra^ and raw/ ran accelerated until 
the gel was formed. During rigid DLCA the slowing down 
of the aggregation was not observed so that gelation at 
this volume fraction was faster than for flexible DLCA. 
The growth rate at different volume fractions can be un- 
derstood by considering the structure of the clusters as 
will be discussed below. 

Figure(l2]) shows N{ra) for ra^] = 100 obtained at dif- 
ferent volume fractions. At the lowest concentration we 
observed a bell shape as is also observed for rigid DLCA 
in the fiocculation regimei^s^i. At higher volume frac- 



3 




10-4 lQ-3 lQ-2 £0-1 IQO 10' 102 IQ3 IQ4 

t 




102 103 10-* 



FIG. 1: (a) rrijn is plotted as a function of time t obtained 
from rigid DLCA (filled symbols) and flexible DLCA (open 
symbols) for e = 0.1 at different (j> as indicated in the figure. 
The solid line has a slope of 1. (b) The polydispersity index 
is plotted as a funtion of t for different (j>- 



tions N(m) followed approximately the power law decay 
expected for the percolation regime: N{m) cx m~^'^ . We 
note that the cut-off function of the size distribution of 
percolating clusters has a Gaussian shapei^i^. More ex- 
tensive simulations on even larger boxes will be needed 
to test whether N{m) obtained from flexible DLCA has 
exactly the same functional form as for rigid DLCA or 
static percolation. 

Aggregation leads to an increase of the number of 
bonds per particle (z). In fig. ([3^) the average number 
of bonds per particle (< z >) is plotted as a function of 
time. As mentioned above, the starting value of < z > 
is larger for larger (p because more randomly distributed 
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FIG. 2: Cluster size distribution is plotted as a function of 
the aggregation number m obtained from fiexible DLCA for 
e — 0.1 at different as indicated in the figure. The solid line 
has a slope of —2.2. 



spheres are within each others interaction range. At low 
volume fractions, < z > increased during rigid DLCA 
until it stagnated at a value close to 2, while at large (j) 
the increase of < z > was small because almost all the 
particles were already part of the final structure at the 
start of the aggregation. During flexible DLCA < z > in- 
creased sharply at first, but when it had reached a value 
around 7 further growth of < z > became very slow. 
Contrary to the case of rigid DLCA, < z > continued to 
increase over the whole duration of the simulation even 
after all spheres were attached to the network indicating 
that restructuring of the system occurred. This process 
was very slow and may be called ageing. For </> > 0.08 
most of the restructuring happened after the gel point 
and continued even when the sol fraction had become 
very small. 



In fig.lISb) are plotted the values of < z > at the end 
of the simulation as a function of the volume fraction. 
< z > has a minimum at « 0.2, the origin of which will 
be discussed below. The distribution of z at long times 
only weakly depended on the volume fraction for < 0.3. 
It was approximately Gaussian with a half width of 4. It 
started at z = 3 and peaked close to < z >. Very few 
particles had z — 12 showing that no crystallisation had 
occurred. At higher volume fractions the distribution 
shifted to larger z. 
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FIG. 3: (a) The average number of neighbors obtained from 
rigid DLCA (filled symbols) and flexible DLCA (open sym- 
bols) for e = 0.1 is plotted as a function of time for different 
(j> as indicated in the figure. Crosses indicate the gel time for 
flexible DLCA. (b) < 2 > at the end of the simulation for 
flexible DLCA with e = 0.1 is plotted as a function of cj>. 



B. Cluster structure 



the value of < z > at the gel point was close to the final 
value for (f) < 0.05, but much lower at large cj), see fig.ijSK). 

Figure (jl^) shows images of clusters formed during flex- 
ible DLCA at very low volume fractions. The configura- 
tion of clusters up to m = 6 was unique contrary to clus- 
ters formed by rigid DLCA. For m ~ 7 one configuration 
was observed in 90% of the cases, but one other config- 
uration was also possible. With increasing aggregation 
number more different configurations were found. Large 
clusters were generally elongated strands that branched 
when m exceeded about 20. The randomly branched 
structures resembled clusters formed by rigid DLCA, ex- 
cept that the strands were much thicker than the single 
particle diameter, see fig. (HJd) . The basic unit of all larger 
clusters formed by flexible DLCA was a tetrahedron and 
for m > 4 particles with less than 3 neighbours were 
not observed. The fact that the tetrahedron is the basic 
structural unit does not imply however that tetrahedra 
were formed first and subsequently aggregated to form 
larger clusters as was suggested ivX. At all stages of the 
aggregation we observed all aggregation numbers with no 
preference for multiples of 4. 
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FIG. 4: (a)Images of clusters with different aggregation num- 
bers formed by flexible DLCA with e = 0.1. (b) Comparison 
of a cluster formed by flexible DLCA (left) and by rigid DLCA 
(right) with aggregation number 700. 



During rigid DLCA the cluster configuration is deter- 
mined by random collisions and the clusters cannot rear- 
range once they are formed. Therefore rigid DLCA yields 
in the flocculation regime clusters with the same average 
radius of gyration for a given aggregation number inde- 
pendent of time. This is not the case for flexible DLCA 
where the particles in the clusters rearrange until a max- 
imum of bonds are formed without breaking any existing 
bond. In this way the density of the clusters increases 
with time. At very low volume fractions the growth rate 
was slower than the restructuring time so that at each 
moment the cluster configuration had reached the steady 
state, while at high volume fractions the gel was formed 
quicker than the time needed to restructure. Therefore 



As mentioned in the introduction, the large scale struc- 
ture of clusters formed by rigid DLCA in the flocculation 
regime is self similar with c?/ = 1.8. In fig. ((5^) a com- 
parison is shown of the dependence of m on Rg between 
clusters formed by flexible DLCA and by rigid DLCA at 
(j) = 0.01. In both cases m = a ■ -Rg ^, but the prefac- 
tor was larger for flexible DLCA (a = 8.5) than for rigid 
DLCA (a = 3.8) because the local structure of clusters 
is denser. FigurelOa) shows the dependence of m on Rg 
for clusters formed at (p — 0.2 close to the gel point. 
The structure of the clusters was almost the same for 
flexible and rigid DLCA because not much restructuring 
had yet occurred before the gel was formed. For larger 
clusters m oc R^^^ as predicted by the percolation the- 
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ory. At intermediate volume fractions the density of the 
clusters increased with time for flexible DLCA and the 
transition between flocculation and percolation shifted to 
larger values of m with decreasing (jj. 
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FIG. 5: The aggregation number is plotted as function of Rg 
for both rigid DLCA (triangles) and flexible DLCA (circles) 
at (a) = O.Of and (b) = 0.2 with e = O.f . The solid lines 
and the dashed line have slopes of —1.8 and —2.5 respectively. 

The structure of the clusters at higher volume frac- 
tions was the same for rigid and flexible DLCA, yet the 
growth rate was faster for flexible DLCA. This can be 
explained by the fact that neighbouring fractal clusters 
have a larger probability to bind when they are flexible 
than when they are rigid. It is likely that in the case of 
rigid clusters rotation would also accelerate the growth 
during the percolation regime. At low volume fractions 
both the collision radius and the hydrodynamic radius 



of clusters formed during flexible DLCA were smaller for 
a given aggregation number. Apparently, both effects 
compensated to give similar growth kinetics for flexible 
and rigid DLCA. However, the percolation regime was 
reached at larger to during flexible DLCA than during 
rigid DLCA. Densification of the clusters during flexible 
DLCA explains why the growth rate at = 0.05 was 
initially faster and at later times slower than for rigid 
DLCA. Initially, the growth was faster due to flexibil- 
ity, but later densification rendered the clusters smaller 
so that the transition to the percolation regime was de- 
layed. 



C. Gel structure 

The gel structure can be characterized by the pair cor- 
relation function {gir)^ or its Fourier transform the static 
structure factor (S[q)). g{r) represents the average num- 
ber concentration of particles at a distance r from any 
given particle and reaches the number concentration of 
the system (C = (pQ/n) at large r. An extensive study 
of gels formed by rigid DLCA was reported ini^iii. It 
was shown that pair correlation functions of gels at low 
volume fractions had distinct features at small r followed 
by a power law decay {g{r) oc A'^f~^'> with dj = 1.8) for 
r > 3. g(r) had a weak minimum at a characteristic value 
fmin close to the correlation length of the concentration 
fluctuations before it reached C. The correlation length 
decreased with increasing volume fraction, and became of 
the order of a few particle diameters for cj) > 0.05. As a 
consequence, the fractal structure did not exist at higher 
volume fractions. 

g{r) showed a delta peak at r = 1 representing the 
contribution of the two nearest neighbours and increased 
continuously to 0.18 at r = 2 starting from a low value 
close to r = I. At low volume fractions the main contri- 
bution to g(r) in this range came from the bound next- 
nearest neighbours. Their contribution stopped at r = 2 
which caused a discontinuous drop of g{r). A small inflec- 
tion of g{r) was observed at r = 3 marking the influence 
of the second shell and at larger r the power law decay 
started. The local structure of the gels was identical for 
small volume fractions {4> < 0.05) but changed at higher 
volume fractions because positional correlations between 
randomly distributed spheres at t = were no longer 
negligible. 

Figure, ([n^) shows pair correlation hmctions of gels ob- 
tained by flexible DLCA at different volume fractions. 
The values oi g(r) within the interaction range were much 
larger than for r > 1 -I- e and are shown separately in 
fig-lEb)- There is a discontinuity in g(r) at r = I -I- e, i.e. 
the maximum distance between a pair of bonded par- 
ticles. This discontinuity is a consequence of the square 
well interaction and is not seen for a continuous potential 
such as the Lennart-Jones potential. For r > f + e, peaks 
can be seen indicating a high degree of local order that 
was absent for rigid DLCA. The maxima at r = 1.7 — 1.8 
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and r = 2 — 2.2 are characteristic for the tetrahedral 
structure with a bond length between 1 and 1.1. If this 
structure is extended in a regular linear fashion it leads 
to the so-called Bernal spiral. Small peaks at larger dis- 
tances indicate that the order persists to some extent but 
for r > 5 the structure became self similar. 
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FIG. 6: (a) The pair correlation function is plotted for differ- 
ent <j) for e = 0.1 as indicated in the figure. The dashed line 
indicate the interaction range 1.1. (b) Zoom of g(;r) between 
1 < r < 1 + e. 



At even larger distances a very weak minimum was 
observed at small volume fractions similar to that found 
for rigid DLCA, but at a given volume fraction it was 
situated at larger r values for flexible DLCA. At a given 
volume fraction, the correlation length of the gel is larger 
for flexible DLCA because the structure is locally denser. 
For the same reason the radius Rc at the crossover be- 
tween flocculation and percolation was larger. It was 
shown inW, for rigid DLCA, that the correlation length 



and Rc are close. 

Bonded nearest neighbours can be situated at any dis- 
tance within the interaction range, i.e. between 1 and 1.1 
in the present case, but figldb) shows that the distance 
was not uniformly distributed. In fact there was a pref- 
erence for distances close to 1 or 1 + e, which means that 
some bonds were compressed and others stretched. The 
local structure was independent of the volume fraction for 
< 0.1, because it is mainly determined by the aggrega- 
tion process of particles that were initially outside each 
others range. At higher volume fractions a significant 
number of particles were in contact before the aggrega- 
tion process started. The restructuring was thus more 
constrained and leading to less local order. The split be- 
tween the peaks at r = 1.7 and 2 became less distinct 
and resembled more closely the split peak observed for 
super cooled liquids. 

The corresponding structure factors are shown in 
fig-® for different volume fractions. At small q a max- 
imum was found at a value [qmax) that increased with 
increasing concentration. For rigid DLCA it was shown 
that the position of maximum is inversely proportional 
to that of the minimum of the pair correlation function 
<lmax oc 2>/rmirr^- Over a narrow g-range S{q) decreased 
with increasing q following a power law, which is expected 
for self similar structures: S{q) cx q~'^' . The data are 
compatible with dj = 1.8, see the solid line in fig.©, 
but for a precise determination of the fractal dimension 
even smaller volume fractions need to be investigated. At 
large g-values S{q) oscillated around unity. 
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FIG. 7: The structure factor is plotted for e at different </> as 
indicated in the figure. The solid lines are guide to the eye 
and dashed line have a slope of 1.8. For clarity we have kept 
only few symbols. 

Figure® compares the structure factors of stable gels 
formed by flexible and rigid DLCA at = 0.05. The 
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maximum has a much larger ampUtude and is situated 
at smaller g-values for gels formed by flexible DLCA and 
the oscillations at large q are stronger. The reason is 
that the local density of the gels is higher and therefore 
the correlation length is larger. At this volume fraction a 
large scale fractal structure can be observed if the gels are 
produced by flexible DLCA, but not if they are produced 
by rigid DLCA. 
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FIG. 8: The structure factor obtained from flexible DLCA 
(solid lines) and rigid DLCA (dashed line) for (j) — 0.05 and 
e = 0.1. 



D. Structural evolution 

The evolution of g{r) and S{q) as a function of time 
during rigid DLCA was shown in^**. With time the max- 
imum of S{q) shifted to lower q and increased in am- 
plitude, while the minimum of g{r) shifted to larger r, 
because the correlation length increased. During rigid 
DLCA the structure is fixed as soon as the bonds are 
formed, but this is not the case for flexible DLCA. We 
illustrate this point by showing the evolution of the struc- 
ture a,t (j) — 0.08 during flexible DLCA. At this volume 
fraction the local structure is not yet influenced by cor- 
relation between the randomly distributed spheres. 

g{r) is plotted in fig® at different times. At the start 
of the aggregation, g{r) was constant within the bond 
range, i.e. 1 < r < 1.1, and equal to C. With time the 
number of bonds increased, but g{r) remained constant 
in the range even when the gel was formed. At longer 
times, however, g(r) became larger at r = 1 and r — 
1.1, implying that in order to maximize the number of 
bonds under constraint, many bonds needed to be either 
stretched or compressed during restructuring of the gel. 



The peaks at r = 1.7 and r = 2 appeared already before 
the gel point indicating that some tetrahedral structure 
had been formed without bond stretching or compression. 
The characteristic minimum of g(r) at rmin developed 
slowly and moved to longer distances. 
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FIG. 9: (a) g{r) is plotted as a funtion of distance for = 0.08 
and e = 0.1 between l<r<l + eat different times as 
indicated in the flgure. (b) g{r) is plotted for r > 1 + e for 
the same system at different times as indicated in the flgure. 

Figure ([TU|) shows the corresponding evolution of the 
structure factor. The shift of the maximum to lower q 
and the increase of its amplitude are caused by an in- 
crease of the correlation length. The amplitude of the 
oscillations at high g- values increased due to local densifi- 
cation. At low volume fractions the pair correlation func- 
tion and the structure factor changed very little after the 
gel point because the aggregation was sufficiently slow 
so that the clusters had time to restructure before they 
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percolated. At high volume fractions gels were formed 
very rapidly and the structure at the gel point resem- 
bled that of gels formed by rigid DLCA. In this case, the 
distinguishing features of the flexible DLCA gels devel- 
oped after the gel point. At long times the structural 
changes became increasingly slow, but they did not stop 
completely during the simulation time. 
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FIG. 10: The S{q) is plotted for <j} = 0.08 and e = 0.1 for 
different times as indicated in the figure. 



E. Dynamics 

For rigid DLCA the mean squared displacement 
(MSD) of the spheres in the gels is, of course, zero. How- 
ever, in flexible DLCA gels the spheres have a significant 
mobility, which can be characterised by measuring the 
MSD as a function of time. It is important to distinguish 
displacements due to restructuring from displacements 
due to flexibility. The displacement due to restructuring 
can be probed by measuring the average total MSD from 
the start of the aggregation process. The total MSD as 
a function of time increases rapidly at first followed by a 
very weak increase when all particles are part of the gel. 

One can characterise the fiexibility of the gels by mea- 
suring the MSD over a period of time during which re- 
structuring is negligible, i.e. for gels that have aged for a 
long time. Alternatively, one can stop further bond for- 
mation and thus further restructuring. This is, of course, 
easier to do in simulations than in experiments on real 
systems. Figure([TT^) shows the MSD of spheres in gels at 
the latest simulation time, when all particles were part of 
the percolating network. Initially, free diffusion was ob- 
served over very short distances. Then the average MSD 
slowed down and finally stagnated at a value < >= 6^. 



S"^ increased with decreasing concentration, but we could 
not properly observe at low volume fractions because 
the times scales involved were too large. 
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FIG. 11; (a) The MSD of spheres in aged gel is plotted as a 
funtion of time for different as indicated in the figure, (b) 
The corresponding distribution of at the longest time. 

If the displacement of the particles is caused by Brow- 
nian motion, the distribution of (or the self part of the 
van Hove correlation function P(r^)) can be described by 
a Gaussian function. Figure pil D) shows P(r^) when the 
average was close to . It is clear that the displace- 
ment was highly heterogeneous as might be expected 
from the heterogeneous structure of the gels. The non- 
Gaussian character can be expressed by the parameter 
a —< r"' > /(< >)^ — 5/3 which is zero for a Gaussian 
function. We found awl for all volume fractions. 

The self part of the intermediate scattering function 
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was found to fully decay at scattering wave vectors q ^ 
1/(5, while it showed a plateau at long delay times for 
smaller g at a value that increased with decreasing q. 
This is, of course, a direct consequence of the restricted 
MSD of the particles, see below. 



F. Effect of the interaction range. 

For rigid DLCA the effect of increasing the interaction 
range is twofold. In the first place, as mentioned above, 
the collision radius of the spheres is increased by a factor 
1 + e. As a consequence the clusters are larger and the 
aggregation is faster. In the second place, the concentra- 
tion of bonds between the randomly distributed spheres 
at the start of the aggregation is larger. The latter effect 
becomes important when the average distance between 
nearest neighbours (A) is smaller than 1 + e (A < 1.1 
for (j> > 0.2). The percolation threshold of randomly 
distributed spheres decreases with increasing e was dis- 
cussed in^. However, as long as A is much larger than 
1 + e, i.e. at small volume fractions, the effect of the in- 
teraction range on the structure of the clusters and the 
gels is small and the fractal dimension remains the same. 

For flexible DLCA the same effects are present, but in 
addition the system has a larger degree of freedom to in- 
crease the number of bonds by restructuring if e is larger. 
As a consequence the value of < z > of the gels increases 
with increasing e. For e = 0.5 we found < z >w 9.5 
after long times at low volume fractions. A qualitatively 
different local structure is formed when e is larger than 
0.41 because the tetrahedral configuration is no longer 
the basic unit. Consequently all spheres are bound to all 
other spheres for pentamers and hcxamcrs. Figure (fT^ 
shows a comparison of the structure factor a,t <j) — 0.05 
for e = 0.1 and e = 0.5. The maximum is situated at 
larger q for e = 0.5 and has a smaller amplitude, indi- 
cating that the correlation length is smaller if e is larger. 
The second maximum is situated at smaller q indicating 
that the local structure (the strand thickness) is larger 
and the oscillations at high q have a lower amplitude in- 
dicating less order. The comparison shows that when the 
interaction range is larger then the gels are more homo- 
geneous with thicker strands and less local order. The 
effect of the interaction range on the structure is small 
for small e and we expect that using a smaller interaction 
range than 0.1 will have little influence on the structure. 



IV. DISCUSSION 

Experimentally, structures similar to those found in 
the present simulations of flexible DLCA have been re- 
ported for mixtures of hard spheres and polymers in 
which an effective attraction between the spheres is in- 
duced by depletion of the polymers^^. If the concentra- 
tion of polymers is high the attraction may become so 
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FIG. 12: The S{q) for e = 0.1 (solid line) and e = 0.5 (dashed 
line) at (j)^ 0.05. 



strong that the aggregation is irreversible on the time 
scale of observation. In this case no significant evolu- 
tion of the structure is observed after gelation. We note 
that in experiments matters are sometimes complicated 
by electrostatic interaction. 

Recently, experiments on aggregating emulsion 
droplets were reported where the bonds were claimed 
to be truly irreversible, but bound droplets were free 
to move along the interface^. The structure factor 
obtained from the experiments was compared to results 
of flexible DLA simulations^. Many of the features 
that distinguish the structure of clusters formed by rigid 
DLCA from those formed by flexible DLCA have also 
been observed in simulations of rigid and flexible DLA 
with zero interaction range (note that the expressions 
classic and slippery were used ini^ to refer to rigid 
and flexible, respectively). Also in DLA, flexible bonds 
caused the formation of thicker strands with a local 
tetrahedral configuration. As for DLCA, the fractal 
dimension of DLA clusters was the same for flexible 
and rigid bonds. The distribution of z was similar 
for flexible DLA and DLCA, and the structure factor 
showed prominent oscillations at high (/-values in both 
cases. 

In fig. (fT5|) the structure factor obtained from flexible 
DLCA (e — 0.1) is compared to that obtained from flex- 
ible DLA (e = 0) and experiments. The experimental 
structure factor was obtained by dividing the scatter- 
ing data of the system with the particle form factor of 
the droplets. The latter was obtained experimentally be- 
fore the aggregation had started and was assumed to 
be not influenced by interaction. As mentioned above, 
the aggregation process is more realistically described 



10 



by DLCA, which results in a smaller fractal dimension 
{df = 1.8) than obtained by DLA {df = 2.5). Nev- 
ertheless, the results from the two types of simulations 
are remarkably close at large g- values, indicating that al- 
most the same local structure is formed. This is perhaps 
not surprising if one considers that local restructuring oc- 
curred in both cases until the maximum number of bonds 
was formed under almost the same constraints. The sim- 
ulations reproduced the position of the peak at q = 7.5 
seen in the experiments, but the amplitude was some- 
what smaller (notice that the experimental data have 
been shifted upward for the comparison with the simula- 
tion results in figure 4 ofi^). The upturn at low g- values 
started at larger q in the simulations than in the exper- 
iments. However, one should be careful when drawing 
conclusions on the basis of the deviation at low q- values, 
because the experimental results are very sensitive to the 
exact shape of the particle form factor in this g-rangc. 




10° 10' 

q 

FIG. 13: Comparison of S{q) between flexible DLCA (solid 
line), slippery DLA^^ (dotted line) and experiment^ (circles). 

As far as we are aware, flexible DLCA with strictly ir- 
reversible bonds has not been simulated before, but the 
present results may be compared with molecular dynam- 
ics simulations of reversible DLCA in the limit of very 
strong attraction. Fofh et al.'^^ studied the structure and 
dynamics of spheres with a square well interaction with 
width e = 0.005 and two different well depths: —2 and 
—20 kT. The interaction strengths may also be expressed 
in terms of the second virial coefficient^: B2 = —44 and 
B2 = —2.9 • 10^, in units of the particle volume. At 
both interaction strengths the systems are far below the 
binodal of the liquid-liquid phase separation. Crystalli- 
sation was avoided by using a bidisperse distribution of 
spheres with slightly different sizes. At smaller inter- 
action strength the system slowly coarsened with time. 



but at B2 = —2.9 ■ 10^ the bonds may be considered 
irreversible on the time scale of the simulation. This sit- 
uation is thus comparable to flexible DLCA. 

The features of the gels obtained by Fofh et al. were 
very close to the ones reported here. Unfortunately, the 
authors did not show pair correlation functions, which 
would allow a more detailed comparison of the local 
structure. But the images of the clusters and the static 
structure factors of the gels resemble closely those shown 
here. The distribution of the bond coordination num- 
ber was almost the same as the one found in this study, 
but < z > was a bit smaller which can be explained 
by the narrower interaction range. Foffi et al. also ob- 
served a minimum of < 2: > at volume fractions close to 
25%, which they speculated to be related to the critical 
point of the liquid-liquid phase separation. We believe 
that the minimum is caused by a combination of oppos- 
ing effects. At very low volume fractions the clusters 
can restructure while they grow, but with increasing vol- 
ume fraction the growth rate increases and much of the 
restructuring has to occur after the percolation thresh- 
old. This increases the constraints on the restructuring 
and thus lowers < z > with increasing (p. At high vol- 
ume fractions, however, the number of bonds increases 
due to crowding even for randomly distributed spheres. 
This leads to an increase of < z > with increasing (j). The 
latter effect is more important when e is larger, which ex- 
plains the stronger increase of < z > found in the present 
study. 

Foffi et al. argued that the gel was formed by a phase 
separation process interrupted by attractive glass tran- 
sition. The same idea has also been put forward by 
other a^^'^^i'^^'^^i'^°''^-^i'^^ . It is clear, that gelation occurs 
only for interaction strengths where the equilibrium state 
would be phase separated and therefore one might indeed 
call the gel formation an interrupted phase separation. 
However, the suggested origin of the interruption is de- 
batable. The concept of attractive glass was introduced 
to account for the slowing down of the dynamics at high 
concentrations with increasing attraction. Even if attrac- 
tive glass formation would be a useful concept for hard 
spheres with a narrow interaction range at high volume 
fractions, it is not clear what it means in the context of 
the arrested state at low concentrations. Many particles 
are situated at the surface of the network strands and 
cannot be in a glassy state. These surface particles can 
escape from their neighbours by breaking all the bonds. 
This process inevitably leads to coarsening. The arrest 
is only inferred if the observation time is shorter than 
the time needed for surface particles to escape. We have 
done preliminary simulations that showed that the rate 
of coarsening decreased with decreasing B2 following a 
power law and is simply too slow at B2 = —2.9 • 10^ to 
observe on the time scale probed by Foffi et al.. 

The results presented here may also be compared with 
earlier Brownian dynamics simulations of colloidal gels 
formed by monodisperse spheres by d'Arjuzon et al4^. 
In those simulations the spheres interacted with a con- 
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tinuous potential that was steeply repulsive at contact 
and attractive over a range of 0.1. The shape of the po- 
tential was chosen to be close to that of depletion inter- 
action. Above a certain interaction strength the spheres 
crystallised^^, but when the interaction was very strong 
crystallisation was no longer observed for the duration of 
the simulation. A detailed study was done of the dynam- 
ics of this system at a single volume fraction (j) = 0.3 and a 
minimum interaction potential of —8kT. For this interac- 
tion strength the life time of the bonds near the potential 
minimum is very long compared to the duration of the 
simulation. However, stretched bonds were still rapidly 
broken and reformed, because the potential goes to zero 
at the outer limit of the interaction range. In spite of 
the different shape of the potential and the reversibility 
of stretched bonds many features of this simulation are 
similar to those reported here for strictly irreversible ag- 
gregation with the same bond range. It was also found 
that < z > increased rapidly at first followed by a very 
slow increase. A similar distribution of the bond coordi- 
nation number was found at long times and g{r) showed 
similar features notably the two peaks indicating tetra- 
hedral structure. The MSD of spheres in the gel and 
also the distribution of displacements at (/> = 0.3 were 
similar to those reported here. The authors argued that 
displacements larger than one particle diameter were due 
to particles breaking of the network. It is clear from the 
present simulations, however, that even in irreversibly 
bound gels large scale mobility is possible due to coop- 
erative motion. It is probable that at lower volume frac- 
tions the maximum displacement is determined by the 
correlation length of the gels. We note that in a recent 
simulation study of more open gels formed by restricting 
the bond angles and thus avoiding the formation of dense 
strands, it was also noted that the MSD displacement can 
be quite large due to cooperative motion^. 

D'Arjuzon et al. studied the self-intermediate scatter- 
ing function in detail and found a fast decay towards a 
plateau followed by a slow decay to zero. The value of 
the plateau was shown to be directly related to S. The re- 
laxation time of the slow decay increased with increasing 
waiting time and became very slow. Preliminary calcula- 
tions showed similar behaviour for the square well system 
studied here. The similarities between the two systems 
indicate that the reversibility of stretched bonds is not 
an essential feature. We believe that in both systems the 
aging is caused by the formation of more (strong) bonds 
under the constraint that existing (strong) bonds cannot 
break at least for the duration of the simulation. 

Lodge and Heyes^ studied spheres interacting with a 
Lennart-Jones potential with varying interaction range 
and depth using Brownian dynamics simulations. A sim- 
ilar shape of the pair correlation function was found if the 
range was narrow and the interaction strong. However, 
they also observed a peak at a smaller distance indicat- 
ing that some crystallisation had occurred. The authors 
studied the phase separation kinetics up to t = 40 in 



the time units used here. The kinetics was very slow 
in the case of strong attraction and may thus also be 
called aging. In fact the time dependence of the peak 
position shown in fig. pl}]) is close to that obtained by 
Lodge and Heycs using a Lennart-Jones potential with 
a narrow interaction range. This indicates again that 
when the interaction is strong, the breaking of bonds is 
not important for the aging process at least in the early 
stages. 



V. CONCLUSION 

Irreversible DLCA with flexible bonds causes locally 
densification, but on large length scales the structure is 
the same as for DLCA with rigid bonds. The fractal di- 
mension of clusters formed in very dilute systems is 1.8 
in both cases. Locally the systems have a tetrahedral 
structure and show a certain degree of order that is in- 
dependent of the volume fraction for < 0.1. At higher 
volume fraction the order is less distinct and similar to 
that of super cooled liquids. 

The system tries to maximize the number of bonds un- 
der the constraint of no bond breaking. The increase of 
the bond coordination number is fast at first, but it be- 
comes progressively slower. The slow restructuring (ag- 
ing) of the gels continues for very long times. At low 
volume fractions most of the local densification occurs 
while the clusters are formed, but at high volume frac- 
tions {(f> > 0.1) the restructuring occurs mainly after gela- 
tion. 

At very low volume fractions the growth rate of the 
clusters during flexible DLCA is the same as for rigid 
DLCA and can be understood in terms of Smolechowski's 
kinetic equations. At high volume fractions the growth 
rate is somewhat faster, because flexibility increases the 
collision rate. 

Gels formed at low volume fractions show a large de- 
gree of flexibility. The average MSD stagnates at a value 
that increases with decreasing volume fraction. However, 
the MSD of the particles is highly heterogeneous reflect- 
ing the fractal structure of the gels. 

Many features that appeared during irreversible flexi- 
ble DLCA closely resemble simulations and experiments 
reported in the literature on spheres with a strong but fi- 
nite attractive interaction. For the interpretation of these 
studies it is essential to consider whether the effect of 
bond breaking is significant during the time of observa- 
tion. 
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